# module utility

# export cumtrapz, rms, polyder, polyval, polyfit, detrend,
#         loadtxt, readNGA

# using DSP

    function HighPassFilter{T<:Real}(f::T,fs::T;order::Integer=4)
        responsetype = Highpass(f; fs=fs)
        prototype = Butterworth(order)
        digitalfilter(responsetype, prototype)
    end

    function LowPassFilter{T<:Real}(f::T,fs::T;order::Integer=15)
        responsetype = Lowpass(f; fs=fs)
        prototype = FIRWindow(hanning(order+1))
        digitalfilter(responsetype, prototype)
    end

    function BandPassFilter{T<:Real}(f1::T,f2::T,fs::T;order::Integer=4)
        responsetype = Bandpass(f1,f2; fs=fs)
        prototype = Butterworth(order)
        digitalfilter(responsetype, prototype)
    end

    function trapz{T<:Real}(x::AbstractArray{T,1}, y::AbstractArray{T,1}, init::T=0.0)

        len::Int64 = length(y)
        r = init

        for i in 2:len
            r += (x[i]-x[i-1])*(y[i]+y[i-1])/2.0
        end

        return r
    end

    function trapz{T<:Real}(y::AbstractArray{T,1},dx::T=1.0, init::T=0.0)

        len::Int64 = length(y)
        r = init

        for i in 2:len
            r += dx*(y[i]+y[i-1])/2.0
        end

        return res
    end

    function cumtrapz{T<:Real}(x::AbstractArray{T,1}, y::AbstractArray{T,1}, init::T=0.0)

        len::Int64 = length(y)
        res::AbstractArray{T,1} = zeros(len)
        res[1] = init
        r = init

        for i in 2:len
            r += (x[i]-x[i-1])*(y[i]+y[i-1])/2.0
            res[i] = r
        end

        return res
    end

    function cumtrapz{T<:Real}(y::AbstractArray{T,1},dx::T=1.0, init::T=0.0)

        len::Int64 = length(y)
        res::AbstractArray{T,1} = zeros(len)
        res[1] = init
        r = init

        for i in 2:len
            r += dx*(y[i]+y[i-1])/2.0
            res[i] = r
        end

        return res
    end

    function rms{T<:Real}(a::AbstractArray{T,1})
        return sqrt(1.0/length(a)*sum(a.*a))
    end

    function polyder{T<:Real,S<:Integer}(p::AbstractArray{T,1},m::S=1)
        n = length(p)-1
        y = p[1:end-1].*(n:-1:1)
        if m == 0
            res = p
        else
            res = polyder(y, m-1)
        end
        return res
    end

    function polyval{T<:Real}(p::AbstractArray{T,1},x::T)
        res = 0.0
        for i in 1:length(p)
            res = x*res+p[i]
        end
        return res
    end

    function polyval{T<:Real}(p::AbstractArray{T,1},x::AbstractArray{T,1})
        res = zeros(length(x))
        for i in 1:length(p)
            res = x.*res+p[i]
        end
        return res
    end

    function polyfit{T<:Real,S<:Integer}(a::AbstractArray{T,1},t::AbstractArray{T,1},oh::S=1,ol::S=0)
        n = length(a)
        A = zeros(n,oh-ol+1)

        for i = 1:oh-ol+1
            A[:,i] = t.^(oh-i+1)
        end

        [A\a;zeros(ol)]
    end

    function leastsq{T<:Real}(A::AbstractArray{T,2},b::AbstractArray{T,1})
        B = deepcopy(A)
        return Base.LinAlg.LAPACK.gelsd!(B,b)
    end

    function detrend{T<:Real,S<:Integer}(a::AbstractArray{T,1},oh::S=1,ol::S=0)
        n = length(a)
        t = linspace(0.0,(n-1)*0.01,n)
        p = polyfit(a,t,oh,ol)
        a - polyval(p,t)
    end

    function detrend!{T<:Real,S<:Integer}(a::AbstractArray{T,1},oh::S=1,ol::S=0)
        n = length(a)
        t = linspace(0.0,(n-1)*0.01,n)
        p = polyfit(a,t,oh,ol)
        a = a - polyval(p,t)
        nothing
    end

    function basecorrd{T<:Real,S<:Integer}(acc::AbstractArray{T,1},disp::AbstractArray{T,1},t::AbstractArray{T,1},oh::S=4,ol::S=2)
        pd = polyfit(disp,t,oh,ol)
        pa = polyder(pd,2)
        acc - polyval(pa,t)
    end

    function basecorrv{T<:Real,S<:Integer}(acc::AbstractArray{T,1},vel::AbstractArray{T,1},t::AbstractArray{T,1},oh::S=3,ol::S=1)
        pd = polyfit(vel,t,oh,ol)
        pa = polyder(pd,1)
        acc - polyval(pa,t)
    end

    function basecorrd{T<:Real,S<:Integer}(acc::AbstractArray{T,1},disp::AbstractArray{T,1},dt::T=1.0,oh::S=4,ol::S=2)
        n = length(acc)
        t = linspace(0.0,dt*n-dt,n)
        pd = polyfit(disp,t,oh,ol)
        pa = polyder(pd,2)
        acc - polyval(pa,t)
    end

    function basecorrv{T<:Real,S<:Integer}(acc::AbstractArray{T,1},vel::AbstractArray{T,1},dt::T=1.0,oh::S=3,ol::S=1)
        n = length(acc)
        t = linspace(0.0,dt*n-dt,n)
        pd = polyfit(vel,t,oh,ol)
        pa = polyder(pd,1)
        acc - polyval(pa,t)
    end

    function beginwithzero{T<:Real}(acc::AbstractArray{T,1})

        Na = length(acc)
        dt = 1.0
        t = linspace(0.0,dt*Na-dt,Na)

        a0 = acc[1]
        aslope = (0.0-a0)/(t[end]-t[1])

        acc = acc .- aslope.*(t .- t[1]) .- a0

        return acc
    end

    function endwithzero{T<:Real}(acc::AbstractArray{T,1},vend::T,dend::T)

        f(i::Int64,t::Float64)   = t^(i+1+1)
        fp(i::Int64,t::Float64)  = i+1>-1 ? (i+1+1)*t^(i+1): 0.0
        fpp(i::Int64,t::Float64) = i+1>0  ? (i+1+1)*(i+1)*t^(i+1-1): 0.0

        Na = length(acc)
        dt = 1.0
        t = linspace(0.0,dt*Na-dt,Na)
        tc = t[end] # controlled time
        b = zeros(3) # right-hand-side vector
        A = zeros(3,3) # coeff matrix

        # setup b, A
        b[1] = 0.0 - dend
        b[2] = 0.0 - vend
        b[3] = 0.0 - acc[end]

        for i = 1:3
            A[1,i] = f(i,tc)
            A[2,i] = fp(i,tc)
            A[3,i] = fpp(i,tc)
        end

        # solve linear system to get c[i]
        c = A\b
        # println((A*c.-b))

        # initial correction time series
        corr = zeros(Na)

        # calculate corr using solved c[i]
        for k = 1:Na
            tk = eqs.t[k]
            for i = 1:3
                corr[k] += c[i]*fpp(i,tk)
            end
        end

        # apply correction
        return acc .+ corr
    end

    function loadtxt(fname::String;skiprows::Integer=0)
        readdlm(fname,Float64;skipstart=skiprows)
    end

    function savetxt{T<:Real}(fname::String,A::AbstractArray{T};delim::String=" ")
        writedlm(fname,A,delim)
    end

    function readNGA(fname::String)
        f = open(fname,"r")
        title = readline(f)
        infor = split(readline(f),",")
        unit = readline(f)
        npts, dt = split(readline(f))[1:2]
        # npts, dt = parse(Int,npts),parse(Float64,dt)
        npts, dt = parseint(npts),parsefloat(dt)
        a_str = replace(readall(f),"\n"," ")
        a_spl = split(a_str)
        # a = [parse(Float64,ai) for ai in a_spl]
        a = Float64[parsefloat(ai) for ai in a_spl]
        close(f)
        return a, npts, dt
    end

    function zero_cross{T<:Real}(a::AbstractArray{T,1})

        # 零交点个数
        nzc = 0
        for i in 1:(length(a)-1)
            if a[i]*a[i+1]<0.0
                nzc += 1
            end
        end
        return nzc

    end

    function peak{T<:Real}(a::AbstractArray{T,1})
        return maxabs(a)
    end

    function peakRatio{T<:Real}(a::AbstractArray{T,1},a0::AbstractArray{T,1})
        return peak(a)/peak(a0)
    end

    function newmarkBeta{T<:Real}(GA::AbstractArray{T,1},DT::T,P::T,ZETA::T;BETA::T=0.25,GAMMA::T=0.5)

        N    = length(GA)
        # TIME = zeros(T,N)
        RD   = zeros(T,N)
        RV   = zeros(T,N)
        RA   = zeros(T,N)

        W = 6.2831854/P

        B1 = 1.0/(BETA*DT*DT)
        B2 = 1.0/(BETA*DT)
        B3 = 1.0/(BETA*2.0)-1.0
        B4 = GAMMA/(BETA*DT)
        B5 = GAMMA/BETA-1.0
        B6 = 0.5*DT*(GAMMA/BETA-2.0)
        B7 = DT*(1.0-GAMMA)
        B8 = DT*GAMMA

        K = W*W
        C = 2.0*ZETA*W

        KEFF = K+B1+B4*C
        KINV = 1.0/KEFF

        RL = [0.0,0.0,0.0]
        RC = [0.0,0.0,0.0]

        RP = [0.0,0.0,0.0] # peak responses
        RT = [0.0,0.0,0.0] # coresponding time

        for I = 1:N

            GC = GA[I]

            FEFF  = GC+(B1*RL[1]+B2*RL[2]+B3*RL[3])+C*(B4*RL[1]+B5*RL[2]+B6*RL[3])
            RC[1] = FEFF*KINV
            RC[2] = B4*(RC[1]-RL[1])-B5*RL[2]-B6*RL[3]
            RC[3] = GC-K*RC[1]-C*RC[2]

            RL[1] = RC[1]
            RL[2] = RC[2]
            RL[3] = RC[3]

            # TIME[I] = I*DTf

            RD[I] = -RC[1]
            RV[I] = -RC[2]
            RA[I] = -RC[3]+GC # absolute acc

        #     if PEAKONLY
        #         if abs(RA[I])>abs(RP[1])
        #             RP[1] = RA[I]
        #             RT[1] = I
        #         end
        #         if abs(RV[I])>abs(RP[2])
        #             RP[2] = RV[I]
        #             RT[2] = I
        #         end
        #         if abs(RD[I])>abs(RP[3])
        #             RP[3] = RD[I]
        #             RT[3] = I
        #         end
        #     end
        end

        # if PEAKONLY
        #     return RP,RT
        # else
        #     return RA,RV,RD
        # end
        return RA,RV,RD
    end

    function RungeKutta{T<:Real}(GA::AbstractArray{T,1},DT::T,P::T,ZETA::T)
        # SUBROUTINE RUNGE_KUTTA(T,RD,RV,RA,GA,DT,ZETA,P,N)
        # !
        # !  SUBROUTINE FOR COMPUTATION OF RESPONSE OF SDOF UNDER EARTHQUAKE
        # !  GA -----  Ground Acc
        # !  N  -----  length of GA
        # !  DT  -----  Original TIME INTERVAL
        # !  ZETA  ----  Damping Ratio

        N = length(GA)
        RD   = zeros(T,N)
        RV   = zeros(T,N)
        RA   = zeros(T,N)

        W=2.0*pi/P

        TW = 2.0*ZETA*W
        W2 = W*W

        H=DT
        HD2=H/2.0
        HD6=H/6.0
        # !
        # !  INITIATION
        # !
        RL = [0.0,0.0,0.0]
        RC = [0.0,0.0,0.0]
        G  = [0.0,0.0]
        GL = 0.0
        # !
        # !  COMPUTATION OF RESPONSE
        # !
        for I=1:N

            G[1] = GL
            G[2] = GA[I]

            A1 = G[1]
            A2 = (G[2]+G[1])*0.5
            A3 = A2
            A4 = G[2]

            X1 = RL[1]
            Y1 = RL[2]
            F1 =-A1+RL[3]

            X2 = RL[1]+Y1*HD2
            Y2 = RL[2]+F1*HD2
            F2 =-A2-TW*Y2-W2*X2

            X3 = RL[1]+Y2*HD2
            Y3 = RL[2]+F2*HD2
            F3 =-A3-TW*Y3-W2*X3

            X4 = RL[1]+Y3*H
            Y4 = RL[2]+F3*H
            F4 =-A4-TW*Y4-W2*X4

            RC[1]=RL[1]+HD6*(Y1+2.0*Y2+2.0*Y3+Y4)
            RC[2]=RL[2]+HD6*(F1+2.0*F2+2.0*F3+F4)

            RC[3]=-TW*RC[2]-W2*RC[1]

            RL[1] = RC[1]
            RL[2] = RC[2]
            RL[3] = RC[3]

            GL = GA[I]
            RD[I] = RC[1]
            RV[I] = RC[2]
            RA[I] = RC[3]
        end
        return RA,RV,RD
    end


    function stepIntegrate{T<:Real}(acc::AbstractArray{T,1},dt::T,P::T,zeta::T)

        r = int(50*dt/P)
        if r>1
            dt = dt/r
            acc = sinc_interp(acc,r)
        end

        ra,rv,rd = newmarkBeta(acc,dt,P,zeta)

        if r>1
            ra = resample(ra,r)
            rv = resample(rv,r)
            rd = resample(rd,r)
        end

        return ra,rv,rd
    end

    function duhamelConv{T<:Real}(acc::AbstractArray{T,1},dt::T,P::T,zeta::T;RA_ONLY::Bool=false)

        t = linspace(0.0,dt*(length(acc)-1),length(acc))
        omega_n = 2.*pi/P
        omega_d = omega_n*sqrt(1.-zeta^2)
        ont = omega_n*t
        odt = omega_d*t
        sinodt = sin(odt)
        cosodt = cos(odt)
        expontz = exp(-ont*zeta)
        if RA_ONLY
            h_dd = expontz.*(-2.*omega_n*zeta*cosodt+(omega_n^2*zeta^2/omega_d-omega_d)*sinodt)
            u_dd = conv(-acc,h_dd*dt)
            ra = u_dd[1:length(acc)]
            return ra
        else
            h = expontz.*sinodt/omega_d
            h_d = expontz.*(cosodt-omega_n*zeta*sinodt/omega_d)
            h_dd = expontz.*(-2.*omega_n*zeta*cosodt+(omega_n^2*zeta^2/omega_d-omega_d)*sinodt)

            u    = conv(-acc,h*dt)
            u_d  = conv(-acc,h_d*dt)
            u_dd = conv(-acc,h_dd*dt)

            rd =    u[1:length(acc)]
            rv =  u_d[1:length(acc)]
            ra = u_dd[1:length(acc)]

            return ra,rv,rd
        end
    end

    function freqResponse{T<:Real}(acc::AbstractArray{T,1},dt::T,P::T,zeta::T;RA_ONLY::Bool=false)

        N = length(acc)
        Nfft = nextpow2(nextpow2(N)+1)
        af = fft([acc;zeros(Nfft-N)])

        w0 = 2.0*pi/P
        f = fftfreqs(Nfft,1.0/dt)
        w = 2.0*pi.*f

        if RA_ONLY
            Ha = zeros(Complex128,Nfft)

            for i = 1:Nfft
                Ha[i] = (w0^2+2.0*zeta*w0*w[i]*1.0im)/(w0^2-w[i]^2+2.0*zeta*w0*w[i]*1.0im)
            end

            ra = real(ifft(af.*Ha))[1:N]

            return ra
        else
            Hd = zeros(Complex128,Nfft)
            Hv = zeros(Complex128,Nfft)
            Ha = zeros(Complex128,Nfft)

            for i = 1:Nfft
                Hd[i] = -1.0/(w0^2-w[i]^2+2.*zeta*w0*w[i]*1.0im)
                Hv[i] = -w[i]*1.0im/(w0^2-w[i]^2+2.0*zeta*w0*w[i]*1.0im)
                Ha[i] = (w0^2+2.0*zeta*w0*w[i]*1.0im)/(w0^2-w[i]^2+2.0*zeta*w0*w[i]*1.0im)
            end

            rd = real(ifft(af.*Hd))[1:N]
            rv = real(ifft(af.*Hv))[1:N]
            ra = real(ifft(af.*Ha))[1:N]

            return ra,rv,rd
        end
    end

    function findZeroCross{T<:Real}(a::AbstractArray{T,1})
        nzc = 0
        zc = Int64[]
        for i = 1:length(a)
            if a[i]*a[i-1]<=0.0
                nzc += 1
                abs(a[i])<=abs(a[i-1]) ? push!(zc,i) : push!(zc,i-1)
            end
        end
        return zc,nzc
    end

    function peakResponseSDOF{T<:Real}(acc::AbstractArray{T,1},dt::T,P::T,zeta::T;method::String="mixed",interp::Function=lin_interp,np::Integer=20)

        r = 1
        if method=="newmark" || method=="N"
            r = int(ceil(np*dt/P))
            if r>1
                dt = dt/r
                acc = interp(acc,r)
            end
            ra,rv,rd = newmarkBeta(acc,dt,P,zeta)
        elseif method=="duhamel" || method=="D"
            r = int(ceil(np*dt/P))
            if r>1
                dt = dt/r
                acc = interp(acc,r)
            end
            ra,rv,rd = duhamelConv(acc,dt,P,zeta)
        elseif method=="freq" || method=="F"
            ra,rv,rd = freqResponse(acc,dt,P,zeta)
        elseif method=="mixed" || method=="M"
            if P<np*dt
                ra,rv,rd = freqResponse(acc,dt,P,zeta)
            else
                ra,rv,rd = newmarkBeta(acc,dt,P,zeta)
            end
        else
            warn("Invalid method, using newmark-beta method.")
            r = int(ceil(np*dt/P))
            if r>1
                dt = dt/r
                acc = sinc_interp(acc,r)
            end
            ra,rv,rd = newmarkBeta(acc,dt,P,zeta)
        end

        if r > 1
            ra,rv,rd = downsample(ra,r),downsample(rv,r),downsample(rd,r)
        end
        ida,idv,idd = indmax(abs(ra)),indmax(abs(rv)),indmax(abs(rd))
        pra,prv,prd = ra[ida],rv[idv],rd[idd]

        return pra,prv,prd,ida,idv,idd
    end

    function AccResponseSDOF{T<:Real}(acc::AbstractArray{T,1},dt::T,P::T,zeta::T;method::String="newmark",INTP::Bool=true,np::Integer=10,interp::Function=lin_interp)

        if method=="newmark" || method=="N"
            if INTP
                r = int(ceil(np*dt/P))
                if r>1
                    dt = dt/r
                    acc = interp(acc,r)
                end
            end
            ra, = newmarkBeta(acc,dt,P,zeta)
            if INTP
                ra = downsample(ra,r)
            end
        elseif method=="runge_kutta" || method=="R"
            if INTP
                r = int(ceil(np*dt/P))
                if r>1
                    dt = dt/r
                    acc = interp(acc,r)
                end
            end
            ra, = RungeKutta(acc,dt,P,zeta)
            if INTP
                ra = downsample(ra,r)
            end
        elseif method=="duhamel" || method=="D"
            ra = duhamelConv(acc,dt,P,zeta;RA_ONLY=true)
        elseif method=="freq" || method=="F"
            ra = freqResponse(acc,dt,P,zeta;RA_ONLY=true)
        elseif method=="mixed" || method=="M"
            if P<np*dt
                ra = freqResponse(acc,dt,P,zeta;RA_ONLY=true)
            else
                ra, = newmarkBeta(acc,dt,P,zeta)
            end
        else
            warn("Invalid method, using newmark-beta method.")
            r = int(ceil(np*dt/P))
            if r>1
                dt = dt/r
                acc = interp(acc,r)
            end
            ra, = newmarkBeta(acc,dt,P,zeta)
            ra = resample(ra,r)
        end

        return ra
    end

    function findpeaks{T<:Real}(a::AbstractArray{T,1})
        
        Na = length(a)
        peaks_pos = zeros(T,Na)
        pinds_pos = ones(Integer,Na)
        peaks_neg = zeros(T,Na)
        pinds_neg = ones(Integer,Na)

        dap = a[2] - a[1]
        np_pos = 0
        np_neg = 0
        for i = 2:Na-1
            da = a[i+1] - a[i]
            if dap >= 0.0 && da <= 0.0 && a[i]>0
                np_pos += 1
                peaks_pos[np_pos] = a[i]
                pinds_pos[np_pos] = i
            elseif dap <= 0.0 && da >= 0.0 && a[i]<0
                np_neg += 1
                peaks_neg[np_neg] = a[i]
                pinds_neg[np_neg] = i
            end
        end

        if np_pos == 0 && np_neg == 0
            warn("no peak found...")
            return nothing,nothing
        elseif np_pos == 0
            return peaks_neg[1:np_neg],pinds_neg[1:np_neg]
        elseif np_neg == 0
            return peaks_pos[1:np_pos],pinds_pos[1:np_pos]
        else
            return peaks_pos[1:np_pos],pinds_pos[1:np_pos],peaks_neg[1:np_neg],pinds_neg[1:np_neg]
        end
    end

    function getEnvelop{T<:Real}(a::AbstractArray{T,1},n::Integer)
        
        peaks_pos,pinds_pos,peaks_neg,pinds_neg = findpeaks(a)
        
        for i = 1:n
            peaks_pos, pinds_pos_tmp = findpeaks(peaks_pos)
            peaks_neg, pinds_neg_tmp = findpeaks(peaks_neg)
            pinds_pos = pinds_pos[pinds_pos_tmp]
            pinds_neg = pinds_neg[pinds_neg_tmp]
        end

        Na = length(a)

        evlp_ind_pos = [1;pinds_pos;Na]
        evlp_pos = [a[1];pinds_pos;a[end]]
        evlp_ind_neg = [1;pinds_neg;Na]
        evlp_neg = [a[1];pinds_neg;a[end]]

        return evlp_pos,evlp_ind_pos,evlp_neg,evlp_ind_neg
    end

    function moveAverage{T<:Real}(SPA::AbstractArray{T,1})
        ASPA = similar(SPA)
        ASPA[1] = SPA[1]#*0.5+0.25*SPA[2]+0.25*SPA[3]
        ASPA[end] = SPA[end]#*0.5+0.25*SPA[end-1]+0.25*SPA[end-2]
        ASPA[2:end-1] = T[ 0.25*SPA[i-1] + 0.5*SPA[i] + 0.25*SPA[i+1] for i=2:length(SPA)-1 ]
        return ASPA
    end

    function moveAverage{T<:Real}(SPA::AbstractArray{T,1},n::Integer)

        ASPA = deepcopy(SPA)
        for i = 1:n
            ASPA = moveAverage(ASPA)
        end
        return ASPA
    end

    function codeSpectrum{T<:Real}(Tg::T,SPT::AbstractArray{T,1},zeta::T=0.05;amax::T=2.25)

        gamma = 0.9+(0.05-zeta)/(0.3+6.0*zeta)
        eta1  = 0.02+(0.05-zeta)/(4.0+32.0*zeta)
        eta2  = 1.0+(0.05-zeta)/(0.08+1.6*zeta)

        SP0 = 1.0/amax
        SPA = ones(length(SPT)).*eta2

        idx1 = searchsorted(SPT,0.1).stop
        idx2 = searchsorted(SPT,Tg).start
        idx3 = searchsorted(SPT,5.0*Tg).stop

        SPA[1:idx1] = SP0+10.0*(eta2-SP0)*SPT[1:idx1]
        SPA[idx2:idx3] = eta2*(Tg./SPT[idx2:idx3]).^gamma
        SPA[idx3+1:end] = eta2*0.2^gamma-eta1*(SPT[idx3+1:end]-5.0*Tg)

        SPA = amax*SPA
        SPV = SPA.*SPT/(2.*pi)
        SPD = SPV.*SPT/(2.*pi)

        return SPA,SPV,SPD

    end

    function codeSpectrum{T<:Real}(Tg::T,SPT::AbstractArray{T,1},Zeta::AbstractArray{T,1};amax::T=2.25)

        NT = length(SPT)
        NZ = length(Zeta)
        SPA = zeros(T,NT,NZ)
        SPV = zeros(T,NT,NZ)
        SPD = zeros(T,NT,NZ)

        for k = 1:NZ
            SPA[:,k],SPV[:,k],SPD[:,k] = codeSpectrum(Tg,SPT,Zeta[k];amax=amax)
        end

        return SPA,SPV,SPD

    end

    function responseSpectrum{T<:Real}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},dt::T=0.02,zeta::T=0.05;method::String="mixed",ifpseudo::Bool=false,ifabs::Bool=true)

        SPA = similar(SPT)
        SPV = similar(SPT)
        SPD = similar(SPT)

        INDSPA = ones(Integer,length(SPT)) # index of peak response (acc)

        for i = 1:length(SPT)
            P = SPT[i]

            if P <= 0.0
                SPA[i] = peak(acc)
                SPV[i] = 0.0
                SPD[i] = 0.0
                continue
            end

            ra,rv,rd,ida,idv,idd = peakResponseSDOF(acc,dt,P,zeta;method=method)
            SPD[i] = ifabs? abs(rd) : rd
            if ifpseudo
                SPV[i] = SPD[i]*(2.*pi)/P
                SPA[i] = SPV[i]*(2.*pi)/P
                INDSPA[i] = idd
            else
                SPV[i] = ifabs? abs(rv) : rv
                SPA[i] = ifabs? abs(ra) : ra
                INDSPA[i] = ida
            end
        end
        return SPA,SPV,SPD,INDSPA
    end

    function responseSpectrum{T<:Real}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},dt::T,Zeta::AbstractArray{T,1};method::String="mixed",ifpseudo::Bool=false,ifabs::Bool=true)

        NT = length(SPT)
        NZ = length(Zeta)
        SPA = zeros(T,NT,NZ)
        SPV = zeros(T,NT,NZ)
        SPD = zeros(T,NT,NZ)

        INDSPA = ones(Integer,NT,NZ) # index of peak response (acc)

        for k = 1:NZ
            SPA[:,k],SPV[:,k],SPD[:,k],INDSPA[:,k] = responseSpectrum(acc,SPT,dt,Zeta[k];method=method,ifpseudo=false,ifabs=false)
        end
        return SPA,SPV,SPD,INDSPA
    end

    function fftConj!{T<:Complex}(A::AbstractArray{T,1})
        Nfft = length(A)
        if iseven(Nfft)
            A[Nfft/2+2:end] = reverse(conj(A[2:Nfft/2]))
        else
            A[(Nfft-1)/2+2:end] = reverse(conj(A[2:(Nfft-1)/2+1]))
        end
        return nothing
    end

    function responseSimple{T<:Real}(Omega::T,omega::T,phi::T,zeta::T,time::AbstractArray{T,1})
        res = similar(time)
        for (i,t) = enumerate(time)
            res[i] = exp(-omega*zeta*t)*((4*omega*Omega^3*zeta-4*omega*Omega^3*zeta^3)*exp(omega*zeta*t)*sin(Omega*t+phi)+((2*Omega^4-2*omega^2*Omega^2)*zeta^2-2*Omega^4+2*omega^2*Omega^2)*exp(omega*zeta*t)*cos(Omega*t+phi)+sqrt(4*omega^2-4*omega^2*zeta^2)*(4*cos(phi)*omega*Omega^2*zeta^3+2*sin(phi)*Omega^3*zeta^2+(cos(phi)*omega^3-3*cos(phi)*omega*Omega^2)*zeta-sin(phi)*Omega^3+sin(phi)*omega^2*Omega)*sin(sqrt(4*omega^2-4*omega^2*zeta^2)*t/2.0)+(8*cos(phi)*omega^2*Omega^2*zeta^4+4*sin(phi)*omega*Omega^3*zeta^3+(2*cos(phi)*omega^4-10*cos(phi)*omega^2*Omega^2)*zeta^2-4*sin(phi)*omega*Omega^3*zeta+2*cos(phi)*omega^2*Omega^2-2*cos(phi)*omega^4)*cos(sqrt(4*omega^2-4*omega^2*zeta^2)*t/2.0))/(8*omega^2*Omega^2*zeta^4+(2*Omega^4-12*omega^2*Omega^2+2*omega^4)*zeta^2-2*Omega^4+4*omega^2*Omega^2-2*omega^4)
            res[i] = -res[i] + cos(Omega*t+phi)
        end
        return res
    end

    function responseSimple{T<:Real}(Omega::T,omega::T,phi::T,zeta::T,t::T)
        return cos(Omega*t+phi) - exp(-omega*zeta*t)*((4*omega*Omega^3*zeta-4*omega*Omega^3*zeta^3)*exp(omega*zeta*t)*sin(Omega*t+phi)+((2*Omega^4-2*omega^2*Omega^2)*zeta^2-2*Omega^4+2*omega^2*Omega^2)*exp(omega*zeta*t)*cos(Omega*t+phi)+sqrt(4*omega^2-4*omega^2*zeta^2)*(4*cos(phi)*omega*Omega^2*zeta^3+2*sin(phi)*Omega^3*zeta^2+(cos(phi)*omega^3-3*cos(phi)*omega*Omega^2)*zeta-sin(phi)*Omega^3+sin(phi)*omega^2*Omega)*sin(sqrt(4*omega^2-4*omega^2*zeta^2)*t/2.0)+(8*cos(phi)*omega^2*Omega^2*zeta^4+4*sin(phi)*omega*Omega^3*zeta^3+(2*cos(phi)*omega^4-10*cos(phi)*omega^2*Omega^2)*zeta^2-4*sin(phi)*omega*Omega^3*zeta+2*cos(phi)*omega^2*Omega^2-2*cos(phi)*omega^4)*cos(sqrt(4*omega^2-4*omega^2*zeta^2)*t/2.0))/(8*omega^2*Omega^2*zeta^4+(2*Omega^4-12*omega^2*Omega^2+2*omega^4)*zeta^2-2*Omega^4+4*omega^2*Omega^2-2*omega^4)
    end

    function responseSimple{T<:Real}(omega::T,phi::T,zeta::T,t::T)
        return cos(Omega*t+phi) + exp(-omega*zeta*t)*(sqrt(4*omega^2-4*omega^2*zeta^2)*(exp(omega*zeta*t)*sin(omega*t+phi)+(-2*cos(phi)*zeta-sin(phi))*cos(sqrt(4*omega^2-4*omega^2*zeta^2)*t/2.0))+(4*cos(phi)*omega*zeta^2+2*sin(phi)*omega*zeta-2*cos(phi)*omega)*sin(sqrt(4*omega^2-4*omega^2*zeta^2)*t/2.0))/(zeta*sqrt(4*omega^2-4*omega^2*zeta^2))/2.0
    end

    function Envelope{T<:Real}(t::AbstractArray{T,1},Tb::T,Tc::T)

        Td = t[end]

        a = -log(0.1)/(Td-Tc)
        N = length(t)
        E = ones(T,N)

        for i = 1:N
            ti = t[i]
            if ti < Tb
                E[i] = (ti/Tb)^2
            elseif ti > Tc
                E[i] = e^(-a*(ti-Tc))
            end
        end

        return E

    end

    function Envelope{T<:Real}(t::AbstractArray{T,1})

        Td = t[end]
        Tb = 0.2*Td
        Tc = 0.5*Td

        return Envelope(t,Tb,Tc)

    end

    function RandPhi{T<:Real}(evlp::AbstractArray{T,1})

        N = length(evlp)
        Nfft = nextpow2(N)

        E = [evlp;zeros(T,Nfft-N)]

        dt = 1.0
        tcum = linspace(0.0,dt*Nfft-dt,Nfft)
        Ecum = cumtrapz(tcum,E)
        Ecum = Ecum./Ecum[end]
        tcum = tcum./tcum[end]

        # plot(tcum,Ecum)

        Rn = rand(int(Nfft/2))
        dphi = lin_interp(Ecum,tcum,Rn)*2.0*pi
        dphi = reverse(dphi)

        phi = zeros(T,int(Nfft/2+1))

        for i = 2:Nfft/2+1
            phi[i] = phi[i-1] + dphi[i-1]
        end

        mod(phi,2.0*pi),dphi

    end

    function findfirstpeak{T<:Real}(a::AbstractArray{T,1})
        da = 0.0
        dap = 0.0
        ida = 1
        while da*dap>=0.0
            dap = da
            da = a[ida+1] - a[ida]
            ida += 1
        end
        return ida
    end

    function findlastpeak{T<:Real}(a::AbstractArray{T,1})
        da = 0.0
        dap = 0.0
        ida = length(a)
        while da*dap>=0.0
            dap = da
            da = a[ida-1] - a[ida]
            ida -= 1
        end
        return ida
    end

    function ArtificialWave{T<:Real}(t::AbstractArray{T,1},evlp::AbstractArray{T,1},SPT::AbstractArray{T,1},SPAT::AbstractArray{T,1},zeta::T=0.05;ifdiffphase::Bool=false,method::String="mixed",tol::T=0.05,maxiter::Integer=15,printout::Bool=true)

        N = length(t)
        dt = t[2]-t[1]
        Nfft = nextpow2(N)

        if ifdiffphase
            phi, = RandPhi(evlp)
        else
            phi = rand(int(Nfft/2+1))*2.0*pi
        end

        freqs = fftfreqs(Nfft,1.0/dt)
        PF = 1.0./freqs[1:Nfft/2]
        Pidx = find(x->x>=SPT[1]&&x<=SPT[end],PF)
        PC = PF[Pidx]

        SPAF = lin_interp(SPT,SPAT,PC)
        p = 0.9
        Saw = 2.0*zeta*SPAF.^2./(pi*2*pi*freqs[Pidx].*(-2.0*log(-pi*log(p)./(2.0*pi*freqs[Pidx]*t[end]))))

        Ak = 2.0*sqrt(Saw*2*pi*(1.0/dt)*Nfft/2)

        af = zeros(Complex128,Nfft)

        af[Pidx] = Ak.*(cos(phi[Pidx]).+1im.*sin(phi[Pidx]))
        fftConj!(af)

        a = real(ifft(af)[1:N].*evlp)
        oh=4; ol=2
        v = cumtrapz(a,dt)
        a = basecorrv(a,v,dt,max(oh-1,0),max(ol-1,0))
        v = cumtrapz(a,dt)
        d = cumtrapz(v,dt)
        a = basecorrd(a,d,dt,oh,ol)
        a = filt(HighPassFilter(1.0/(SPT[end]),1.0/dt),a)
        v = cumtrapz(a,dt)
        d = cumtrapz(v,dt)
        fpk = findfirstpeak(v)
        a[1:fpk] = beginwithzero(a[1:fpk])
        lpk = findlastpeak(v)
        lpk = findlastpeak(v[1:lpk-1])
        a[lpk:end] = endwithzero(a[lpk:end],v[end],d[end])

        a,SPA,SPV,SPD,INDSPA = SPFit(a,SPT,SPAT,dt,zeta;method=method,tol=tol,maxiter=maxiter,printout=printout)
        return a,SPA,SPV,SPD,INDSPA

    end

    function adjustpeak!{T<:Real}(a::AbstractArray{T,1},peakvalue::T=1.0)
        pka = peak(a)

        if pka>peakvalue
            a[find(x->x>peakvalue,a)] = peakvalue
            a[find(x->x<-peakvalue,a)] = -peakvalue
        end

        if pka<peakvalue
            idmx = indmax(a)
            idmn = indmin(a)
            if a[idmx]>-a[idmn]
                a[idmx] = peakvalue
            else
                a[idmn] = -peakvalue
            end
        end
    end

    function SPFit{T<:Real}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},SPAT::AbstractArray{T,1},dt::T=0.02,zeta::T=0.05;method::String="mixed",tol::T=0.05,maxiter::Integer=15,printout::Bool=true,Nci::Integer=-1,Ncf::Integer=-1)

        NT = length(SPT)
        Na = length(acc)

        if Nci<0
            Nci = findfirst(x->abs(x)>0.1*peak(acc),acc)
            Nci = max(Nci,int(Na*0.1))
        end

        if Ncf<0
            Ncf = findfirst(x->abs(x)>0.1*peak(acc),reverse(acc))
            Ncf = max(Ncf,int(Na*0.1))
        end

        ExpFac = 0.5

        NSPT = [SPT[1]*(1.0-ExpFac);SPT;SPT[end]*(1.0+ExpFac)]
        NSPAT = [SPAT[1]-(SPAT[2]-SPAT[1])/(SPT[2]-SPT[1])*SPT[1]*ExpFac;SPAT;SPAT[end]+(SPAT[end]-SPAT[end-1])/(SPT[end]-SPT[end-1])*SPT[end]*ExpFac]

        SPA,SPV,SPD,INDSPA = responseSpectrum(acc,NSPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)

        R = abs(NSPAT./SPA)
        Rm1 = abs(1.0./R .- 1.0)
        aerror = sqrt(sum(Rm1[2:end-1].^2)/(NT-2))
        merror = maximum(Rm1[2:end-1])

        acc, Nfft, Na = fft_padding(acc;base=2)
        af = fft(acc)
        ai = deepcopy(acc)

        freqs = fftfreqs(Nfft,1.0/dt)
        PF = 1.0./freqs[1:Nfft/2]
        Pidx = find(x->x>=NSPT[1]&&x<=NSPT[end],PF)
        PC = PF[Pidx]

        RF = lin_interp(NSPT,R,PC)

        # oh=4
        # ol=2

        iter = 1
        while (aerror>tol || merror>2.0*tol) && iter<=maxiter-1
            if printout
                println("Iter = $iter, Error = $(round(aerror,3)), MaxErr = $(round(merror,3))")
            end
            # modify the multiplier acorrding to the sign of every freq response at the time of wave's peak response
            j = 1
            for i = 1:NT
                P = SPT[i]
                pkt = dt*INDSPA[i]-dt
                omega = 2.0*pi/P
                if i == 1
                    dP = P-NSPT[i]
                elseif i == NT
                    dP = NSPT[i+1]-P
                else
                    dP = min(P-NSPT[i],SPT[i+1]-P)
                end

                while PC[j]>P-0.5*dP && PC[j]<P+0.5*dP
                    Omega = 2.0*pi/PC[j]
                    phi = angle(af[Pidx[j]])
                    rai = responseSimple(Omega,omega,phi,zeta,pkt)
                    if sign(rai)*sign(SPA[i]) < 0.0
                        RF[j] = 1.0/RF[j]
                    end
                    j += 1
                end
            end
            # adjust the fourier amplitude to fit target spectrum
            af[Pidx] = RF.*af[Pidx]
            # make a conj mirror against the fold frequency
            fftConj!(af)
            # inverse FFT to recover the acc
            ai = real(ifft(af))

            # remove poly trends in the beginning and the end
            # if iter > 4
            #     ai[1:Nci] = detrend(ai[1:Nci],4,0)
            #     ai[1:Nci] = beginwithzero(ai[1:Nci])
            #     ai[end-Ncf+1:end] = detrend(ai[end-Ncf+1:end],4,1)
            # end

            # ai[1:Nci] = filt(HighPassFilter(10.0/SPT[end],1.0/dt),ai[1:Nci])
            # ai[end-Ncf+1:end] = filt(HighPassFilter(10.0/SPT[end],1.0/dt),ai[end-Ncf+1:end])

            # vi = cumtrapz(ai,dt)
            # di = cumtrapz(vi,dt)
            # ai[end-int(N*0.02):end] = endwithzero(ai[end-int(N*0.02):end],vi[end],di[end])

            # vi = cumtrapz(ai,dt)
            # ai = basecorrv(ai,vi,dt,max(oh-1,0),max(ol-1,0))
            # vi = cumtrapz(ai,dt)
            # di = cumtrapz(vi,dt)
            # ai = basecorrd(ai,di,dt,oh,ol)
            # ai = filt(HighPassFilter(1.0/(4.0*SPT[end]),1.0/dt),ai)
            # ai = filt(BandPassFilter(1.0/(2.0*SPT[end]),1.0/(2.01*dt),1.0/dt),ai)

            adjustpeak!(ai)

            SPA,SPV,SPD,INDSPA = responseSpectrum(ai[1:Na],NSPT,dt,zeta;method=method,ifpseudo=false)

            # figure(10000,(9,6))
            # plot(SPT,abs(SPA[2:end-1]))
            # xscale("log")

            R = abs(NSPAT./SPA)
            RF = lin_interp(NSPT,R,PC)

            Rm1 = abs(1.0./R .- 1.0)
            aerror = sqrt(sum(Rm1[2:end-1].^2)/NT)
            merror = maximum(Rm1[2:end-1])
            iter += 1
        end
        if printout
            println("Iter = $iter, Error = $(round(aerror,3)), MaxErr = $(round(merror,3))")
        end

        return ai[1:Na],abs(SPA[2:end-1]),abs(SPV[2:end-1]),abs(SPD[2:end-1]),INDSPA[2:end-1]
    end

    function SPAdjust{T<:Real}(acc::AbstractArray{T,1},P::T,RAT::T,dW::T,dt::T=0.02,zeta::T=0.05;method::String="freq")

        ra,rv,rd = freqResponse(acc,dt,P,zeta)

        I = indmax(abs(ra))
        RA = ra[I]

        W0 = 2.0*pi/P

        Na = length(acc)
        t = linspace(0.0,dt*Na-dt,Na)
        tp = t[I]
        tmp = t .- tp
        ap = sign(RA).*(RAT-abs(RA))*cos(W0.*tmp).*sinc(dW.*tmp)

        apt, Nfft, Na = fft_padding(ap;base=2)
        apf = fft(apt)

        F = fftfreqs(Nfft,1.0/dt)
        W = 2.0*pi.*F
        Ha = (W0^2+2.0*zeta*W0*W*1.0im)./(W0^2-W.^2.+2.0*zeta*W0*W*1.0im)

        return acc .+ real(ifft(apf./Ha))[1:Na]

    end

    function alpha{T<:Real}(f::T,f1::T=1.,f2::T=4.,z1::T=1.25,z2::T=0.25)
        if f<=f1
            tmp = z1*f
        elseif f>=f2
            tmp = z2*f
        else
            tmp = (z1+(z2-z1)*(f-f1)/(f2-f1))*f
        end
        return tmp
    end

    function wavelet_func{T<:Real}(t::AbstractArray{T,1},tm::T,P::T,zeta::T;wtype::Integer=1)

        w = 2.*pi/P
        f = 1./P

        tmp1 = sqrt(1.0-zeta^2)
        gamma = 1.178*(f*tmp1)^-0.93
        deltaT = atan(tmp1/zeta)/(w*tmp1)

        tmp2 = t.-tm+deltaT

        if wtype == 1
            wfunc = cos(w*tmp1*tmp2).*exp(-(tmp2/gamma).^2)
        elseif wtype == 2
            wfunc = cos(w*tmp1*tmp2).*exp(-abs(tmp2)*alpha(f))
        end

        return wfunc
    end

    function SPAdjust{T<:Real}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},SPAT::AbstractArray{T,1},zeta::T=0.05,dt::T=0.02;method::String="mixed",tol::T=0.05,maxiter::Integer=5,printout::Bool=true)

        Na = length(acc)
        t = linspace(0.0,dt*Na-dt,Na)
        SPA,SPV,SPD,INDSPA = responseSpectrum(acc,SPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)
        NT = length(SPT)

        R = abs(SPA./SPAT)
        Rm1 = abs(R .- 1.0)
        aerror = sqrt(sum(Rm1.^2)/NT)
        merror = maximum(Rm1)
        # SR = sortperm(Rm1,by=x->-x)
        # INDA = SR[1:min(NT,50)]
        INDA = 1:NT
        INDSPAP = deepcopy(INDSPA)

        iter = 1
        while (aerror>tol || merror>1.0*tol) && iter<=maxiter-1
            if printout
                println("Iter = $iter, Error = $(round(aerror,3)), MaxErr = $(round(merror,3))")
            end

            acc = spadjust(acc,SPT[INDA],SPA[INDA],INDSPA[INDA],SPAT[INDA],zeta,dt;method=method)

            SPA,SPV,SPD,INDSPA = responseSpectrum(acc,SPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)

            # INDPT = Integer[]
            # for i = 1:NT
            #     if INDSPAP[i] != INDSPA[i]
            #         wf = wavelet_func(t,INDSPA[i]*dt-dt,SPT[i],zeta)
            #         raw = AccResponseSDOF(wf,dt,SPT[i],zeta;method=method)
            #         ra = AccResponseSDOF(acc,dt,SPT[i],zeta;method=method)
            #         dra = (abs(ra[INDSPA[i]])-abs(ra[INDSPAP[i]]))*sign(ra[INDSPA[i]])
            #         acc = acc + 1.0*dra/raw[INDSPA[i]]*wf
            #         append!(INDPT,[i])
            #     end
            # end

            # if length(INDPT)>0
            #     for i = INDPT
            #         ra,rv,rd,ida,idv,idd = peakResponseSDOF(acc,dt,SPT[i],zeta;method=method)
            #         SPA[i] = ra
            #         INDSPA[i] = ida                    
            #     end
            # end

            R = abs(SPA./SPAT)
            Rm1 = abs(R .- 1.0)
            aerror = sqrt(sum(Rm1.^2)/NT)
            merror = maximum(Rm1)
            # SR = sortperm(Rm1,by=x->-x)
            # INDA = SR[1:min(NT,50)]
            INDSPAP = deepcopy(INDSPA)
            iter += 1
        end
        if printout
            println("Iter = $iter, Error = $(round(aerror,3)), MaxErr = $(round(merror,3))")
        end
        return acc,SPA,SPV,SPD,INDSPA
    end

    function spadjust{T<:Real,S<:Integer}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},SPA::AbstractArray{T,1},INDSPA::AbstractArray{S,1},SPAT::AbstractArray{T,1},zeta::T=0.05,dt::T=0.02;method::String="mixed")

        Na = length(acc)
        t = linspace(0.0,dt*Na-dt,Na)
        # SPA,SPV,SPD,INDSPA = responseSpectrum(acc,SPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)
        dR = sign(SPA).*(SPAT.-abs(SPA))./SPAT

        NT = length(SPT)

        WF = AbstractArray{T,1}[wavelet_func(t,INDSPA[j]*dt-dt,SPT[j],zeta) for j = 1:NT]
        A = zeros(T,NT,NT)
        for i = 1:NT
            for j = 1:NT
                wf = WF[j]
                ra = AccResponseSDOF(wf,dt,SPT[i],zeta;method=method)
                A[i,j] = ra[INDSPA[i]]/SPAT[i]
                if i!=j
                    A[i,j] *= 0.618
                end
            end
        end

        b, = leastsq(A,dR)
        da = zeros(Na)
        for i = 1:NT
            da = da .+ b[i]*WF[i]
        end

        # dv = cumtrapz(da,dt)
        # dd = cumtrapz(dv,dt)
        # figure(10000,(12,3))
        # plot(t,dd)

        # SPA,SPV,SPD,INDSPA1 = responseSpectrum(acc,SPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)
        # println(INDSPA1.-INDSPA)

        acc = acc .+ da
        # adjustpeak!(acc)
        # acc = filt(HighPassFilter(1.0/(2.0*SPT[end]),1.0/dt),acc)
        return acc

    end

    function SPAdjust{T<:Real}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},SPAT::AbstractArray{T,2},Zeta::AbstractArray{T,1},dt::T=0.02;method::String="mixed",tol::T=0.03,maxiter::Integer=5,printout::Bool=true)

        NT = length(SPT)
        NZ = length(Zeta)

        SPA,SPV,SPD,INDSPA = responseSpectrum(acc,SPT,dt,Zeta;method=method,ifpseudo=false,ifabs=false)

        R = abs(SPA./SPAT)
        Rm1 = abs(R .- 1.0)
        aerror = sqrt(sum(Rm1[:].^2)/(NT*NZ))
        merror = maximum(Rm1[:])
        # SR = sortperm(Rm1[:],by=x->-x)
        # INDA = SR[1:min(int(NT/2),30)]
        INDA = 1:NT

        iter = 1
        while (aerror>tol || merror>1.0*tol) && iter<=maxiter-1
            if printout
                println("Iter = $iter, Error = $(round(aerror,3)), MaxErr = $(round(merror,3))")
            end

            acc = spadjust(acc,SPT[INDA],SPA[INDA,:],INDSPA[INDA,:],SPAT[INDA,:],Zeta,dt;method=method)

            SPA,SPV,SPD,INDSPA = responseSpectrum(acc,SPT,dt,Zeta;method=method,ifpseudo=false,ifabs=false)

            R = abs(SPA./SPAT)
            Rm1 = abs(R .- 1.0)
            aerror = sqrt(sum(Rm1[:].^2)/(NT*NZ))
            merror = maximum(Rm1[:])
            # SR = sortperm(Rm1,by=x->-x)
            # INDA = SR[1:min(int(NT/2),30)]

            iter += 1
        end
        if printout
            println("Iter = $iter, Error = $(round(aerror,3)), MaxErr = $(round(merror,3))")
        end
        return acc,SPA,SPV,SPD,INDSPA
    end

    function spadjust{T<:Real,S<:Integer}(acc::AbstractArray{T,1},SPT::AbstractArray{T,1},SPA::AbstractArray{T,2},INDSPA::AbstractArray{S,2},SPAT::AbstractArray{T,2},Zeta::AbstractArray{T,1},dt::T=0.02;method::String="mixed")

        # multi-damping-ratio spectra fitting

        Na = length(acc)
        t = linspace(0.0,dt*Na-dt,Na)
        # SPA,SPV,SPD,INDSPA = responseSpectrum(acc,SPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)
        dR = sign(SPA).*(SPAT.-abs(SPA))

        NT = size(SPT,1)
        NZ = length(Zeta)
        NTZ = NT*NZ

        WF = AbstractArray{T,1}[wavelet_func(t,INDSPA[j,k]*dt-dt,SPT[j],Zeta[k]) for j = 1:NT, k = 1:NZ]
        A = zeros(T,NTZ,NTZ)

        for i = 1:NTZ
            for j = 1:NTZ
                wf = WF[j]
                K,I = divrem(i,NT)
                if I == 0
                    I = NT
                    K = K - 1
                end
                ra = AccResponseSDOF(wf,dt,SPT[I],Zeta[K+1];method=method)
                A[i,j] = ra[INDSPA[i]]
                if i!=j
                    A[i,j] *= 0.618
                end
            end
        end

        # for Pi = 1:NT
        #     for Pj = 1:NT
        #         for Zi = 1:NZ
        #             for Zj = 1:NZ
        #                 wf = WF[Pj,Zj]
        #                 ra = AccResponseSDOF(wf,dt,SPT[Pi],Zeta[Zi];method=method)
        #                 A[NZ*(Pi-1)+Zi,NZ*(Pj-1)+Zj] = ra[INDSPA[Pi,Zi]]
        #                 if NZ*(Pi-1)+Zi!=NZ*(Pj-1)+Zj
        #                     A[NZ*(Pi-1)+Zi,NZ*(Pj-1)+Zj] *= 0.618
        #                 end
        #             end
        #         end
        #     end
        # end

        b, = leastsq(A,dR[:])
        da = zeros(Na)
        for i = 1:NT*NZ
            da = da .+ b[i]*WF[i]
        end

        # SPA,SPV,SPD,INDSPA1 = responseSpectrum(acc,SPT,dt,zeta;method=method,ifpseudo=false,ifabs=false)
        # println(INDSPA1.-INDSPA)

        acc = acc .+ da
        # adjustpeak!(acc)
        # acc = filt(HighPassFilter(1.0/(2.0*SPT[end]),1.0/dt),acc)
        return acc

    end

    function Error{T<:Real}(A::AbstractArray{T,1},B::AbstractArray{T,1})
        Rm1 = 1.0 .- A./B

        return sqrt(sum(Rm1.^2)/length(A)),maxabs(Rm1)
    end

    function fft_padding{T<:Real}(a::AbstractArray{T,1};base::Integer=2)
        Na = length(a)
        Nfft = nextpow(base,Na)
        res = zeros(Nfft)
        res[1:Na] = a
        return res, Nfft, Na
    end

    function fftfreqs(Nfft::Integer,fs::Real=1.0)
        fn = 0.5*fs
        df = fs/Nfft
        freqs = zeros(Nfft)
        if iseven(Nfft)
            freqs[2:Nfft/2] = linspace(df,fn-df,int(Nfft/2-1))
            freqs[Nfft/2+1:end] = linspace(-fn,-df,int(Nfft/2))
        else
            l = int((Nfft+1)/2-1)
            freqs[2:(Nfft+1)/2] = linspace(df,df*l,l)
            freqs[(Nfft+1)/2+1:end] = linspace(-df*l,-df,l)
        end

        return freqs
    end

    function sinc_interp_fd{T<:Real}(a::AbstractArray{T,1},r::Integer)

        if r == 1
            return a
        end

        a, Nfft, Na = fft_padding(a)
        af = fft(a)

        Nffti = Nfft*r
        afi = zeros(Complex128,Nffti)

        if iseven(Nffti)
            afi[1:Nfft/2] = af[1:Nfft/2]
            afi[end-(Nfft/2-2):end] = af[end-(Nfft/2-2):end]
            afi[end-(Nfft/2-1)] = 0.5*af[end-(Nfft/2-1)]
        else
            afi[1:(Nfft-1)/2] = af[1:(Nfft-1)/2]
            afi[1:(Nfft-1)/2+1] = 0.5*af[1:(Nfft-1)/2+1]
            afi[end-((Nfft+1)/2-2):end] = af[end-((Nfft+1)/2-2):end]
            afi[end-((Nfft+1)/2-1)] = 0.5*af[end-((Nfft+1)/2-1)]
        end

        ai = r*real(ifft(afi))[1:Na*r-r+1]

        return ai
    end

    function sinc_interp_td{T<:Real}(a::AbstractArray{T,1},r::Integer)

        if r == 1
            return a
        end

        Na = length(a)
        t = linspace(0,Na-1,Na)
        ti = linspace(0,Na-1,Na*r-r+1)

        T[dot(a, sinc((ti[i]-t))) for i in 1:Na*r-r+1]
    end

    function sinc_interp{T<:Real}(a::AbstractArray{T,1},r::Integer;method::String="freq")
        if (method == "time") || (method == "t") || (method == "T")
            return sinc_interp_td(a,r)
        elseif (method == "freq") || (method == "f") || (method == "F")
            return sinc_interp_fd(a,r)
        else
            error("Unknown method: $method.")
        end
    end

    function upsample{T<:Real}(a::AbstractArray{T,1},r::Integer)
        return kron(a, [zeros(r-1);1.0])[r:end]
    end

    function filt_interp{T<:Real}(a::AbstractArray{T,1},r::Integer)

        if r == 1
            return a
        end

        ap = upsample(a,r)

        responsetype = Lowpass(0.5/r; fs=1.0)
        prototype = FIRWindow(hanning(4*r+1))
        filter = digitalfilter(responsetype, prototype)

        ap = r*filt(filter,ap)

        apsft = similar(ap)
        apsft[1:end-(2*r)] = ap[2*r+1:end]
        apsft[end-(2*r-1):end] = 0.0
        return apsft
    end

    function lin_interp{T<:Real}(a::AbstractArray{T,1},r::Integer)

        if r == 1
            return a
        end

        ap = upsample(a,r)

        for i = 1:length(a)-1
            init = a[i]
            slope = (a[i+1] - a[i])/r

            for j = 1:r-1
               ap[r*(i-1)+1+j] = init + slope*j
            end
        end
        return ap
    end

    function lin_interp{T<:Real}(x::AbstractArray{T,1},y::AbstractArray{T,1},xi::AbstractArray{T,1})

        yi = similar(xi)

        for i = 1:length(xi)
            odr = searchsorted(x,xi[i])
            k1, k2 = odr.stop, odr.start
            if k1 == k2
                yi[i] = y[k1]
            elseif k1 == 0
                slope = (y[k1+1]-y[k2+1])/(x[k1+1]-x[k2+1])
                yi[i] = y[k1+1] + slope*(xi[i]-x[k1+1])
            elseif k2 > length(x)
                slope = (y[k1-1]-y[k2-1])/(x[k1-1]-x[k2-1])
                yi[i] = y[k1-1] + slope*(xi[i]-x[k1-1])
            else
                slope = (y[k1]-y[k2])/(x[k1]-x[k2])
                yi[i] = y[k1] + slope*(xi[i]-x[k1])
            end
        end
        return yi
    end

    function resample{T<:Real}(a::AbstractArray{T,1},r::Integer)

        if r == 1
            return a
        end

        Na = length(a)
        Nfft = nextpow(r,Na)

        af = fft([a;zeros(Nfft-Na)])
        # a, Nfft, Na = fft_padding(a;base=r)
        # af = fft(a)

        Nffti = int(Nfft/r)
        afi = zeros(Complex128,Nffti)

        if iseven(Nffti)
            afi[1:Nffti/2] = af[1:Nffti/2]
            afi[end-(Nffti/2-2):end] = af[end-(Nffti/2-2):end]
            afi[end-(Nffti/2-1)] = 1.0*af[end-(Nffti/2-1)]
        else
            afi[1:(Nffti-1)/2] = af[1:(Nffti-1)/2]
            afi[1:(Nffti-1)/2+1] = 1.0*af[1:(Nffti-1)/2+1]
            afi[end-((Nffti+1)/2-2):end] = af[end-((Nffti+1)/2-2):end]
            afi[end-((Nffti+1)/2-1)] = 1.0*af[end-((Nffti+1)/2-1)]
        end

        ai = (1.0/r)*real(ifft(afi))[1:int((Na-1)/r)+1]

        return ai
    end

    function downsample{T<:Real}(a::AbstractArray{T,1},r::Integer)

        if r == 1
            return a
        end

        return a[1:r:end]
    end
# end
